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Filtered Mass Density Function for Design Simulation of 
High Speed Airbreathing Propulsion Systems 

T.G. Drozda, R.M. Sheikhi and P. Givi 
Department of Mechanical and Aerospace Engineering 
State University of New York at Buffalo 
Buffalo, NY 14260-4400 


Abstract 

The objective of this research is to develop and implement new methodology for large eddy simulation 
of (LES) of high-speed reacting turbulent flows. We have just completed two (2) years of Phase 1 of 
this research. This annual report provides a brief and up-to-date summary of our activities during the 

period: September 1, 2000 through August 31, 2001. m 

In the work within the past year, a methodology termed “velocity-scalar filtered density function 
(VSFDF) is developed and implemented for large eddy simulation (LES) of turbulent flows. In this 
methodology the effects of the unresolved subgrid scales (SGS) are taken into account by considering 
the joint probability density function (PDF) of all of the components of the velocity and scalar vectors. 
An exact transport equation is derived for the VSFDF in which the effects of the unresolved SGS 
convection, SGS velocity-scalar source, and SGS scalar-scalar source terms appear in closed form. The 
remaining unclosed terms in this equation are modeled. 

A system of stochastic differential equations (SDEs) which yields statistically equivalent results to the 
modeled VSFDF transport equation is constructed. These SDEs are solved numerically by a Lagrangian 
Monte Carlo procedure. The consistency of the proposed SDEs and the convergence of the Monte 
Carlo solution are assessed by comparison with results obtained by an Eulerian LES procedure in which 
the corresponding transport equations for the first two SGS moments are solved. The unclosed SGS 
convection, SGS velocity-scalar source, and SGS scalar-scalar source in the Eulerian LES are replaced 
by corresponding terms from VSFDF equation. The consistency of the results is then analyzed for a 
case of two dimensional mixing layer. 


Technical Monitor 

Dr. J. Philip Drummond (Hypersonic Propulsion Branch, NASA LaRC, Mail Stop 197, Tel: 757-864-2298) 
is the Technical Monitor of this Grant. 


Personnel 


Dr. Peyman Givi is the PI of this project. One Graduate Research Assistant (RA), Mr. Tomasz G. Drozda, 
is being supported by this Grant. We also acknowledge collaborations with Professor Stephen B. Pope 
(Cornell University) on various aspects of this project. 


1 



1 Introduction 


The probability density function (PDF) approach has proven useful for large eddy simulation (LES) of 
turbulent reacting flows. 7 - 12 ’ 14 ’ 15 ’ 17 ’ 24 ’ 26 ' 30 The formal means of conducting such LES is by consideration 
of the “filtered density function” (FDF) which is essentially the filtered fine-grained PDF of the transport 
quantities. In all previous contributions, the FDF of the “scalar” quantities is considered: Gao and 
O’Brien, 14 Colucci et al , 24 and Reveillon and Vervisch 26 developed a transport equation for the FDF in 
constant density turbulent reacting flows. Jaberi et al. 2S extended the methodology for LES of variable 
density flows by consideration of the “filtered mass density function” (FMDF), which is essentially the mass 
weighted FDF. The fundamental property of the PDF methods is exhibited by the closed form nature of 
the chemical source term appearing in the transport equation governing the FDF (FMDF). This property 
is very important as evidenced in several applications of FDF for LES of a variety of turbulent reacting 
flows. 24, 26-29 However, since the FDF of only the scalar quantities are considered, all of the “hydrodynamic” 
effects are modeled. In all previous LES/FDF simulations, these effects have been modeled via “non- FDF" 
methods. 

The objective of the work conducted this year is to extend the PDF methodology to also include the SGS 
velocity-scalars. This is facilitated by consideration of the joint “velocity-scalar filtered density function 
(VSFDF). With the definition of the VSFDF, the mathematical framework for its implementation in LES 
is established. A transport equation is developed for the VSFDF in which the effects of SGS convection 
are shown to appear in closed form. The unclosed terms in this equation are modeled in a fashion similar 
to those in the Reynolds-averaged simulation (RAS) procedures. A Lagrangian Monte Carlo procedure 
is developed and implemented for numerical simulation of the modeled VFDF transport equation. The 
consistency of this procedure is assessed by comparing the first two moments of the VFDF with those 
obtained by the Eulerian finite difference solutions of the same moments transport equations. 


2 Formulation 

We consider an incompressible (unit density), isothermal, turbulent reacting flow involving N s species. For 
the mathematical description of this flow, the primary transport variables are the velocity vector 
(i = 1, 2, 3), the pressure p(x, t ), and the species’ mass fractions </> Q (x, t) (a = 1, 2, ... , N s ). The equations 
which govern the transport of these variables in space (x t ) and time (f) are 


** = 0 

dx k 

(la) 

du t du k u l dp i d<Jik 

~dt + dx k ~ dxi dx k 

(lb) 

d^a du k 4> a _ dj£ + s * ^ 

dt dx k dx k 

(lc) 

where S a (x,t) denotes the chemical reaction term for species a, and 4> = [<f>u <h, ■ ■ 
scalar array. Assuming a Newtonian flow with Fick’s law of diffusion, we have. 

, ,<f>N a ] denotes the 

/ dui du k \ 

° ik ~ U \dx k dxi) 

(2) 
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( 3 ) 


tq pa 

4 “ ~ r W t 

where v is the kinematic viscosity and is assumed constant, and I ' 1 is the density weighted diffusion 
coefficient of species a. The Schmidt number, Sc, is defined as: T" = 

LES involves the use of the spatial filtering 6,32 


/ +oo 

f{x',t)G{x',x)dx' 

-00 


(4) 


where function G{x',x) is a filter function, ( f{x,t )) is the filtered value of the transport variable f{x,t). 
Applying such a filter to the governing Eq.(la-lc) yields filtered equations 


d{u k ) 

dx k 

d(u i) d{u k )(uj) 


= 0 


dip) d 2 (ui) dr{u k ,ui) 
+ v- 


dt ' dx k dxi dx k dx k dx k 

d{<t>a) | d(u k ){4> a ) = ^ + (SaW<)> 

dt dx k dx k dx k dx k 

where the subgrid scale terms(SGS), r, are in general defined as 

r(aubj) = {aibj) — (a.i)(bj) 

and satisfy following transport equations 

djiuuuj} SKHyii = 

dx k dx k dx k dx k ox k 


(5a) 

(5b) 

(5c) 

( 6 ) 


dt 


, dut duj\ ( dp\ ( dp V 

4wl 3C's^J +T r5ir T r'5;JJ 


drjuk^uh^ 

dx k 


(7) 


The first unclosed term in the square bracket is dissipation and is often labeled as, e ik . The next two terms 
correspond to velocity-pressure scrambling and are labeled in literature as, rij*. 

dT((t> a ,(f>0 ) d(u k )T{4> Q ,<f>p ) va d 2 T(4> a ,4>p) ' djffa) 

— m — + aS r dx t dx t ax k 


r(u k ,4>g) 


d(<P a ) 


dx k 

dripk^g^ </>/?) 

dx k 


2T°t +r(cl>a,Sg)+T((f> 0 ,S a ) ~ 

\ dx k dx k ) J 


( 8 ) 
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where the term in square bracket is the scalar dissipation. 

dr^Uu^g) 

dt 


+ 


+ 

0 ) 

where the first term in square bracket is the velocity-scalar dissipation, next term is scalai-piessure scram- 
bling, and the last two terms are velocity-scalar diffusion and scalar-velocity diffusion terms. All of the 
terms in square brackets are the unclosed terms that need to be modeled. Triple subgrid scale terms, r, 
are generally defined as 

r(aj,6j,Cfc) = (fli&jCfc) — (ai)T(6j, Cfc) — 

— ( bj)r(ai,Ck ) - { Ck)r(ai,bj ) - { ai)(bj)(c k ) (10) 


8(u,)t(u„^) _ /v + rn gM^d- T(ut , Ui) 
dx k V 2 / dx k dx k 


d(4> a ) 

dx k 


T{u k ,(j> a ) 


i/-r° 


r{ui, S a ) - 


djui) 

dx k 


(* + 


pa , ( dui d<t> a \ (. 

r ''(s’bj* r«.i 


+ 


dT(uk,Ui,<f> a ) 


d 2 Ui 

dx k dx k 


)) 


dx k 


A Velocity-Scalar Filtered Density Function 

The “velocity-scalar filtered density function” (VSFDF), denoted from here on as P, is formally defined as 

/-(-00 

6 ( V, t/>; U {x\ y), <j){x', t)) G{x' - x)dx' (11) 

-OO 


3 A I s 

g{y,x/j',U(x l ,y),(f)(x',t)) = ]^[ <5 (Vi — Ui(x, t)) x <5 (xp a ~ <f>a(x,t)) 

i=l Q=1 

(12) 

where S denotes the delta function, V,r/> are the velocity and scalar state vectors respectively. Term g 
represents so called “fine-grained” density 30 and hence Eq.(ll) defines VSFDF as the spatially filtered 
value of the fine-grained density. In other words, VSFDF represents a density in the velocity-scalar space 
of the fluid around x weighted by filter G, at a time t. With the filter properties specified in Eq.(13-17), 
P has all of the properties of the probability density function(PDF). 1 


0 

R^ 

III 

Q 

R^ 

1 

(13) 

G(x) = G(-x) 

(14) 

G(x) > 0, for all x 

(15) 

'-{-OO 

x m G(x)dx = moments exist for rn > 0 

(16) 

-oo 


f+OO 

/ G(x)dx - 1 

J — OO 

(17) 
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Following the procedure of Pope,^ a VSFDF transport equation can be derived by considering a time 
derivative of Eq.(ll). 


dP 9 nr d\ 

di + dT t (VtP) ~ 


+ 


d 2 P dip) dP 


dx k dx k 
d 


+ 


dx k dV k dip a 


[S a (i>)P] + 


dV k 

d 2 


((&"■»>- ! SM- 

dui duj 


dVidVj 

d 2 

dip a dipQ 

d 2 

dVidrpa 


dx k dx k 
p dipg dip 0 
dx k dx k 

dui dipa 


2T 


dx k dx k 


\V,iPjP 

\V,xp)p 

\V,ip)P 


(18) 


where r = r a = v, and consequently Sc a — 1. A quick check of the above equation can be made 
by integrating it according to Eq.(22) to obtain appropriate moment equations. These equations should 
match with Eq.(5a-5c), and Eq.(7-9) 


B Modeled VSFDF Transport Equation 


A closure needs to be provided for the conditional terms on the RHS of VSFDF Eq.(18) to allow for a 
solution. Haworth and Pope 2 proposed that the Generalized Langevin Model (GLM) be employed in this 
task. An extension of their strategy allows us to present a combination of GLM and linear mean square 
estimation(LMSE) 3 ’ 8 ’ 25 model 


d T ( / dp_ 
dV k lAWfc 


d 2 


dlpadlpi 3 


I V,xP 


d (p) \ p 

dx k ) 


d 2 


r^^iv.v 

dx k dx k 


dV x dV } 

d 2 

dVidlPa 


2T 


dx k dx k 
du , dip Q 


dx k dx k 


\V,ip)P 


^ di^djuj) d 2 / 2 p d(ui) &j _ 
dx k dx k dV l dV j dx k dx k dV { 
f) 1 d ^ f 

'dV ^ Gij ^ Vj ~ ^ + 2 C ° € dV k dV k + 


+ o 7 — [C(f> o L0 {ipa {<^ 0 ))/] + 

U'lpa 

r g{0a)gM d 2i l2r ^) g ^) gz I 2 r ^ a) g2/ 

+i dx k dx k dip a dip 0 dx k dx k dV t dip a dx k dx k dip a 


(19) 


Above model offers full consistency with the exact equations, however, as will become more clear in the 
following section, the last line can not be included because it violates scalar field realizability conditions. 
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C Stochastic System 


The most convenient way of solving VSFDF equation is by the Lagrangian Monte Carlo procedure. Eulerian 
Monte Carlo schemes exist, but they are shown to produce excessive artificial diffusion, which greatly 
degrades LES results. 24 Lagrangian Monte Carlo scheme is based on the principle of equivalent systems. 1 
This principle states that two stochastic systems with different instantaneous behaviors may produce 
identical statistics and satisfy the same Fokker-Plank (PDF, FDF) equation. The following is then true: a 
set of SDEs produces one, and only one, Fokker-Plank equation, while a Fokker-Plank equation can result 
from many different sets of SDEs. In this light one may consider a set of stochastic differential equations 
as the most precise way of describing a random process, and its Fokker-Plank analogous to a filter of this 
random process. Lagrangian Monte Carlo scheme is provided by the following general form of a stochastic 
diffusion process 


dX t 

= D x (t)dt + B x {t)dW* (t) 

(20a) 

dUi 

= Df (t)dt + B u (t)dW? (t) + 



+ F% x (t)dWj c (t) + l%+{t)dW+(t) 

(20b) 

dlpa 

= D+(t)dt + B+{t)dW+{t) + 



+ F+f (: t)dW? (t) + F^ (i)dMf (*) 

(20c) 


Above set consists of a standard stochastic drift (D) and diffusion (B) terms. In addition, difiusional cou- 
pling is introduced between variables by adding extra diffusion terms (F). All of the terms (D, B, F) are im- 
plicitly functions of the state variables, that is, Dpt) — D(X.(t),V{t),ip a (t)]t), B(t) = B(X.(t), V{t),ip a (t),t), 
F(t) = F(X(t),V{t),rp a {t)]t). The W terms denote independent Wiener-Levy processes. 16 The corre- 
sponding Fokker-Planck equation is 


df* 

dt 


+ 


+ 


+ 


■ — (D x f*) — (D u /*) — 

dxi { 1 1 ’ dVt [ l 1 1 dip c 

8 MVl + jnsdCm 


(Dir) 

d 


+ 


K*r) 


dxidx , ' ' ' dxidVj ' ' • dx, 

d - (h ux f*\ + - (b uu f*) H — ( b U4> r) 

.. v®« 1 ) + AV ftVx ^ 1 > dVidipj \ l] J ) 

d 


+ 


dVidx, 


dVidVj 

d 


dip 


k(^) + 3 jk(^ r ) + 555 ft 


( 21 ) 


where b tJ = [E(t)r(t)E r (t)] , 16 E is the diffusion matrix and T is the covariance matrix of the Wiener-Levy 
processes, W X ,W U ,& nd W*. 

The above equations implicitly describe a set of moment transport equations. These transport equations 
can be found either by applying ltd formula 16 to the SDEs and applying a filtering operation, or by 
integrating the Fokker-Plank equation. In the later case the integration takes on a form 


<<2(M)) 



Q{V, Ip] X, t)f(V , V»; x, t)dV dip 


( 22 ) 


where Q is a random variable, ( Q ) is its filtered value, and / is its PDF. 
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C.l Fully Consistent System 


Pope 1 and Gicquel 31 have proposed a system of stochastic differential equations (SDEs), based on the 
VSFDF model in Eq.(19). A set of random variables, X, U, <f>, that correspond to the position, velocity and 
scalar value respectively is introduced into a stochastic diffusion process. The unknown coefficients are 
found by comparing the moments of the below set with the exact set of moment equations, i.e. Eq.(5a-9) 


dX l 

dUi 


dlpa 


Uidt + \J v \ (IW^ 

d(j>) , d 2 (ui) r u 

-X b ^2-x — t; b G ij \*j - ( u ])> 

dxi dx k dx k 

+V^^~dW k x + JC^edW? 

dx k 

d 2 (4> a ) 


dt + 


Q dx k dx k 


- C</.a W (^a - ((pa)) + SaW 


dt + 


+v/^r 


d((pa) 

dx k 


dW ] 


x 


(23a) 


(23b) 


(23c) 


The above set of equations contains a model for the unclosed terms in Eq.(18). It should be noted here 
that the first and the last terms in Eq.(23c) violate realizability of the scalar field and are the source of 
the terms appearing on the last line in Eq.(19). In particular, these terms would allow a mixture fraction, 
for example, to extend beyond physical range of 0 to 1. 

Computation of the Fokker- Plank 16 equation, or the PDF(in this case it is also a FDF), for the diffusion 
process in equations above, yields 


% + Jk m) ~ 


d(p) . , v d 2 (ui) ' 


wrm 


- K - Jw;} IQ." W. - <#.» f\-£- IScWW + 


+ 


v x d 2 f , 0 {Uj) o~j uypg] u j_ 

~2 8x k dx k + dxi dxidVj ^ UxUSa dx.'dxidip 


dx k dx k dip a ' dip c 
-d(uj) d 2 f 


dlp a 


d(<p a ) d 2 f 


v-i d(ui) d(uj) d 2 f 


1 


d 2 f 


+ 2 dx k dx k dVrdVj + 2 C ° e dV k dv k 

-d(ui) d(<pa) d 2 f 


t d{4>*) dfo) a 2 / 

"b 2 \ VSa 1 Sq dx k dx k dtpa dipji 


+ 


v\^ll v J 

+ dx k dVidip a 

The moment transport equations resulting from integration Eq.(22) are 

d(u k ) _ Q 

d 2 (ui) 


dx k 

d(ui) d(u k )(uj) 


dt 


dx k 




( y + ^2 - 


dx k dx k 


d(4>a) + d(u k )(ppcP)_ 


( / , _l tq ((h \\ _ £11 


dr(u k ,Ui) 

dx k 

dr{u k ,(p a ) 


dt 


dx k 


(24) 


(25a) 


(25b) 


(25c) 
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The SGS moment equations are 


drjuuUj) + d(u k )r(ui,uj) _ ^ _ r[uk , Uj )gM + 


’ 9^*; 

+ (l>l — 2y/v\V% + V 3 ) 


2 dx k dx k 
d{ui) d{uj) 


dxk dxk 


+ [Gijbr(ujt,Uj) + GjkT{u k ,Ui) + Coe<5ij] - 


dx k 


dr(uk, Uj, Uj) 

dx k 


dx k 


( 26 ) 


dT(<t> a Ap) , d{Uk)T(4>a,<l>0) __ V\ d 2 j(#aM , ■ jgW _ T ( fo) gM + 


dt 


dx k 


2 dx k dx k 


+ (4 - v/i'ii's* - veil's* + V^Sa^j 


dx k 

^ dj^a) d($0} 

dx k dx k 


dx k 


dT{u k ,4> a ,(}>p) 


— [(<?*« + C^urita, 4> 0 )] + T(<f> a , Sp) + r(^, S a ) 


(27) 


dT(Ui,(fia ) d(Uk)T{Ui,<f>g) V\ 5 2 r(uj, 4>a) / 

777 1 o o Q„ \ Ki l 


dt 


dx k 


,d{(f> a ) Mui) 

I— 


2 dx k dx k ' " dx k 

+ («/!- v^-v/^ + n/^T) oT k o^ + 

, X1 , dT{u k ,Ui,<p a ) 

+ [G lk T(u k ,4> a ) - C (j)a uJT(Ui,<p a )\ + T{Ui,b Q ) - ^ 


(28) 


A comparison of above set to the corresponding exact moment equations reveals that the fully consistent 
stochastic system can be achieved if the coefficients in the SDE’s Eq.(30a-30c) are set as follows. 




V\ = V2 = "3 — v a = v S a ,f},... - 2r 



(29) 

A fully consistent stochastic systems becomes: 




dXi 

= Uidt + VtfdW* 



(30a) 

dUi 

~ 

\- 9 ® +2T? { ’? ) + Gij (Vj - (tij))] 

dxi ax k ox k \ 

dt 4“ 



+ V2T d i Ui) dW k x + s/C^edW? 
dx k 



(30b) 

dlpa 

~ 

[2r f ( t Q> C+.uWa-ttaV + SaW) 
dx k oxk J 

dt ~h 



+Vzr d ^iwF 

dx k 



(30c) 
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The Fokker-Planck equation is: 


dj_ 

dt 


+ i^(Vicf) = r 7r-5- + n7>i- + 

dx k dx k dx k dx k dv k 

a(uj)d{u,) aw S*f 

dx k dx k dVidVj dx k dx k 3V l 

+ WJ~ [Qo W (0a — (0a)) /] 

30a 


and the corresponding moment equations obtained by integration Eq.(22) are 


3K) _ g 

dx k 

d{ui) | 3(u fc )(u a ) _ d(p) | r d 2 (^i) dr[u k ,uj) 

dt dx k dxi dx k dx k dx k 

d(4>a) , d(u k )(<p a ) _ d 2 (<p a ) ,.*v dl -{ukM 

-m~ + — tek T d^r k + {a m dx k 


( 31 ) 


(32a) 

(32b) 

(32c) 


= - rh,.^ - .(»„»# + 

dt dx k dx k ox k ox k ox k 


+ [ GikT{u k ,Uj ) + GjkT{u k ,Ui) + CneSij] 


dx k 
3 t(t i k ,Uj,Uj) 

dx k 


(33) 


. a(-t)r(fe,^i =r ^^M _ T{Hl iA a) 2M - T( „. 

dt dx k dx k dx k ox k dx k 

dT{u k ,4>cn4>0) 


— [(C0Q "b C^) wr (0a) 0/j)] + T (0aj Sp) + T { ( fi(h &a) 


dx k 


(34) 


= r ^r{uuM _ _ r ( Ufc ^ a )gM + 

dt dx k dx k dx k ‘ dx k dx k 

^ , , m / r. \ dT(u k ,Ui,(f)a) 

+ [Gi k r(u k , (f) a ) — C0 Q WT(tli, <f> a)] + r ( w z> *S*q) 


(35) 


where the terms in square brackets provide closure to the corresponding unclosed terms in the exact 
Eq.(5a-5c), and Eq.(7-9). 


C.2 VSFDF System for Consistency Assessment 

The stochastic system of SDEs presented in previous section is fully consistent with exact filtered equations 
and capable of simulating convection and reaction terms in a fully closed form. We have showed this 
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mathematically in prior sections, however, it is necessary to also show that numerically, such a system 
indeed produces results comparable with, most commonly used, Eulerian finite difference schemes. This 
task is accomplished by solving a somewhat simpler and realizable system of SDEs 


dX t 
dU i 


dxp a 


= U{dt + V2TdW t x 

m +2r f^L +G ,,(v i -(u i )) 

dxi oxkdxk 

+ Vzr?S^- dW k + y/c^edw? 

ox k 

— -C<f> a U) ('Ipa — (fia)) dt 


dt 4 - 


(36a) 


(36b) 

(36c) 


Above system employs a fully consistent hydrodynamic system and LMSE model for scalar mixing. 
The Fokker-Plank 16 (PDF ) equation for this system is 


df 

dt 


„ d(« i) 9(uj) &f gw ay 

dx k dx k dVidVj dx k dxkdK 

3 1 3^ f 

~ dVi ^ Gi j ^ Vj 2 C ° €Sl3 dVidVj ^ 

+ S~ [C^U (ipa ~{<Pa)) f] 

(jtpn 


and the corresponding moment equations obtained by integration Eq.(22) are 


d(uk) 

dx k 

d(ui ) dju^j Uj) 

dt dxk 

d{pa) dpUk^p q ) 

dt dx k 


0 

d(p) r d 2 (uj) _ dr(ufc,rq) 
dxi dx k dx k dx k 

r d 2 (<ft Q ) _ dT(u k ,<f>g) 
dx k dx k dx k 


(37) 


(38a) 

(38b) 

(38c) 


dr(ui,Uj) 

dt 


+ 


dxk dxkdxk oxk ax k 

^ ^ c 1 dr(uk,Ui,Uj) 

+ Uj) + G jk' r { u k'> u i) S — dxk 


(39) 


dt dx k dx k dx k dx k ox k 


+ 


2r ^) d^)_ _ + c )uT{<f>aM 

dx k dx k 


dx k 


(40) 
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dT{u u (p a ) 

dt 


+ 

+ 


dx k dx k 

[■ G ik T{u k , 0 a ) - C^LOTim, 0 O )] - 


T(ttfc,u») 


5 { 0 o) 


dx k 

dr{u k ^ Uj, <fi a ) 

dx k 


T(u k ,(f>a ) 


djui) 

dx k 


+ 


( 41 ) 


Above set of equations will be solved using a finite difference scheme and the solution will be assessed for 
consistency with a solution from a Lagrangian Monte Carlo solver. 


D Rotta’s Closure and Model Constants 


The GLM and LMSE have provided a stochastic closure for the LES system. However, we have not 
mentioned a proper form of G i} tensor, nor the values of the closure coefficients. Pope 1 ’ 19 has derived an 
appropriate form of Gij and suggested some guidelines for values of the closure coefficients. 

Gij = -ui Q + Sij (42) 

where u = £ , e = C e ^, k = \t (u k , u k ), and A L is the LES filter size. 

The resulting closure for SGS velocity is equivalent to that of Rotta, 19 

- [ t («<. + t ( uj> **) + 2vt te’ fe)] = 

= \Gi k r{u k ,Uj) + Gj k r(u k ,Ui) + Co^ij) = 

= -w (1 + §C 0 ) [r{ui,uj) ~ lkS tJ ] - | eSij (43) 


(Cfc + C^r(K,h) - 2r^^' 


for the SGS scalar we have, 

- [2T«t (^f, gf); = 

and for SGS velocity-scalar we have, 

- [(i/ + r Q )r ( jg;, + T (0a, §£-) + {^~T ~ ) ( T ( Ui » Afe) " T (^ 0 ’ 

= [( G lk r{u k , 0a) - C^uriui, 0 Q )] = -u(\ + § C 0 + C^) r («i, 0a) 
The closure constants as proposed by Pope 19 are 

Co = 2.1 Ct = 1.0 C*= 1.0 


'i 

dx k dx k j J 


(44) 


(45) 


(46) 


E Numerical Procedure 

The solution of the VSFDF equation (Eqn. 18) provides all of the statistical information about velocity- 
scalar field. The most convenient way of solving this equation is by a Lagrangian Monte Carlo procedure. 
In the Lagrangian description, the FDF is represented by an ensemble of statistically identical Monte Carlo 
particles. These particles, in general, do not represent actual fluid particles, but rather serve as carriers 
of information about flow statistics. In a limited range of time scales, however, Monte Carlo particles do 
provide a direct model of fluid particle behavior. 1 


11 



The information carried by the particles is updated through a temporal integration of the SDE s Eq.(36a- 
36c). The simplest way of performing this stochastic simulation is via the Euler approximation 1 

A -+ 1 = X?+(D?)"At+(B x )"m' , 2 (C*)'' < 47i> > 

£/”+' = u;‘ + (Df)" A( + (B a )" (if) 1 '' 2 (('')” 

+ (F$ x )" (At)'/ 2 «*)” (47b) 

* +1 = <* + (bJ)”ai+ (b+)" (ao 1/2 (ci)” 

+ (f*/)” (Al)‘/ 2 (</)" (47c) 

where D, B, and F are all functions of the state variables at a time step n, and C’s are independent 
standardized Gaussian random variables. Higher order numerical schemes, such as Runge-Kutta, are 
also available for stochastic differential equations, 10,21,22 but caution is advised when selecting one for a 
Monte-Carlo simulation of LES. If the diffusion terms, B and F, strongly depend on the state variables, the 
numerical scheme may alter the solution in a way inconsistent with the true nature of undiscretized set of 
SDEs. 16 Numerical scheme must be consistent with Ito and Stratonovich calculus, and Euler approximation 
satisfies this condition. 


The statistical information transported by Monte Carlo particles and evolving according to the SDE’s 
is evaluated by considering a “finite area”, for 2D, and “finite volume”, for 3D, centered at a LES grid 
point. The area/volume is characterized by a length A#. Statistics are computed based on an ensemble 
of approximately N E particles residing inside it at a given time. This type of “ensemble” approach is 
necessary as the probability of finding a single, or multiple for that matter, Monte Carlo partick (s) at a 
given LES grid point is zero. 18 Fig.l illustrates this concept. An alternative to the ensemble averaging is a 
method proposed by Pope. 1 This method uses least-squares continuous cubic splines with continuous first 
and second derivatives. It must be used when the number of Monte Carlo particles is too small to allow for 
smooth resolution of the first and second derivatives from the stochastic solver. Because the statistics are 
computed based on ensemble of Monte Carlo particles statistical errors are introduced into the simulation. 
Ideally, for reliable Eulerian statistics and minimum numerical dispersion, it is desired to minimize the 
size of ensemble domain while maximizing the number of Monte Carlo particles in it. Consequently the 
statistics would tend to their exact LES value. 


(ui)E 


r{ui,Uj) E 


N e 




(n) 


nGAf 


Ne — » oc 
Ae— ►O 


> («»> 


£( 

nGA 




>• r ( Ui,Uj ) 

Ne^oo 

A £ -yO 


(48) 


(49) 


where subscript E indicates a statistical value obtained from the ensemble centered around the LES grid 
point. This procedure is termed “point estimator”, because it ignores the ensemble average variations of 
the statistical value within the ensemble domain and uses the ensemble averaged value at the LES grid 
point only. 

Solution of the Lagrangian stochastic system, as it is presented in Eq.(36a-36c), is not possible because 
SDE’s contain filtered terms that can not be computed by the stochastic solver. To remedy this problem 
a finite difference LES code is developed and coupled with the stochastic solver. This combination yields 
a finite difference-stochastic hybrid LES solver. 
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Figure 1: Concept of ensemble averaging. Shown are three different ensemble domains: 

1(A e = 0.5A ,Ne « 10), 2{A e = A, N E « 40), 3(A fi = 2.0A,1V£; « 160). Black squares indicate LES grid 
points, while small circles are Monte Carlo particles. 


The filtered pressure field, velocity-scalar field and turbulent kinetic energy field appearing in Eqn. 36a- 
36c, are needed by the stochastic solver to advance the system in time. These quantities are computed via 
the finite-difference LES solver by a method based on the “compact parameter” finite difference scheme 
of Carpenter. This is a variant of the Mac Cormack scheme in which fourth-order compact differences 
are used to approximate the spatial derivatives, and second-order symmetric predictor-corrector sequence 
is employed for time discretization. All of the finite difference operations are conducted on a fixed and 
equally spaced grid points. The transfer of information from the LES grid points to the Lagrangian 
particles is accomplished via interpolation. A second-order interpolation scheme is used. The transfer 
of information from the Lagrangian particles to the finite difference solver is accomplished via ensemble 
averaging described earlier in this section. This procedure is shown schematically in Fig. 2. 

As can be seen in Table 2 a certain “redundancy” results from using a finite difference-stochastic hybrid 
LES solver. This is actually very useful as it allows for constant monitoring of the accuracy of the simulated 
results. This “redundancy” will also be used to show that the Lagrangian Monte Carlo solver is consistent 
with its finite difference counterpart. 
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Figure 2: Numerical Procedure Flow Chart for VSFDF Hybrid Solver 
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Table 1: VSFDF Solution Procedure. VSFDF1 refers to actual production model, while VSFDF1 Consis- 
tency refers to the code used to demonstrate consistency with the finite difference LES solver. 


3 Results 

A Flow Simulated 

Consistency simulations of the Lagrangian Monte Carlo solver were conducted for a 2D temporally de- 
veloping mixing layer. A 2D simulation is sufficient to establish consistency and is chosen to preserve 














computational resources. 

The temporal mixing layer consists of two parallel streams traveling in opposite directions with the same 
speeds. Hyperbolic tangent functions are used to initialize both stream-wise velocity and scalar distribu- 
tions profiles, while a form of exponential function is used to initialize the cross-stream velocity profile^ 
The simulations are conducted for 0 < x < L, and -L/2 <y< L/2, where x,y denote stream-wise and 
cross-stream directions, and the length L is specified such that L = 2nN r /a us , where N r is the number 
of successive vortex pairings and a us is the wavelength of the most unstable linear mode corresponding to 
the mean stream-wise velocity profile at initial time. 4 The formation of large scale structures is facilitated 
by introducing small harmonic, phase-shifted, disturbances containing subharmonics of the most unstable 
mode into the stream-wise and cross-stream velocity profiles. This results in formation of two successive 
vortex pairings. For more detailed description of this “forcing” mechanism reader is referred to Ref. 4 
The flow variables are normalized by the half initial vorticity thickness, L r = and half the velocity 
difference across the layer, U T = ^ (<Mo = and (m) is the Reynolds averaged value of filtered 

stream-wise velocity ) 


B Numerical Specifications 


All of the simulations are performed on the equally spaced finite difference grid, such that Ax = Ay = A. 
The LES resolution used was 32x41 grid points. The choice for resolution was guided by the proper ratio 
of the resolved to total Reynolds stresses. 23 In general for LES, this ratio should be approximately 80%. 
The Reynolds number used for the simulations was, Re = = 50. A constant LES top-hat filter of size 

A i = 2A centered around a grid point is used such that, 


G{x' - x) 


G(x[ - xi) 


No _ 

= G(x' - Xi) 
i=l 



, \ x i - %i\< ^2 

,W-*i l>Tf 


(50) 

(51) 


where Np is the number of dimensions. 

The Monte Carlo particles are initially distributed according to the uniform PDF scaled to fit the domain. 
All of these particles have equal “weight” in the solution. It is required that the particle distribution 
be approximately uniform to minimize the ensemble average bias and interpolating errors. The periodic 
boundary condition in stream-wise direction ensures that the particles leaving the domain are introduce 
at the opposite boundary with the same compositional values. The cross-stream direction has a free-slip 
boundary condition such that leaving particles are “mirror-reflected” back into the domain. The particle 
density is determined by counting the Monte Carlo particles residing inside the ensemble domain of size 
A e x Ae(xAe). The effects of average particle density (N E ) and ensemble domain size (A E ) on the LEb 
solution will be assessed to ensure consistency and statistical accuracy of the VSFDF simulations. 

The results will be analyzed both instantaneously, by visually comparing the finite difference and stochastic 
LES results at a specific time, and statistically, by Reynolds averaging the results at a specific time. In 
the case of a temporally developing shear layer the Reynolds average statistics are constructed by spatially 
averaging over the homogeneous direction and are denoted by an overbar. 
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C Consistency and Convergence Assessments 

To show that the VSFDF formulation in consistent and that the Monte-Carlo simulation converges we 
will compare results obtained via VSFDF and LES-FD. This approach is valid because the accuracy and 
reliability of the finite-difference LES procedures is well established and consequently will provide a good 
means of assessing the quality of the Monte-Carlo simulation of the VSFDF . The comparison will consist 
of statistical and instantaneous values of all of the variables computed by the hybrid solver. No attempt 
is made to determine the appropriate values of the model coefficients, rather the values suggested in the 
literature are adopted (See section D). 

Fig. 3 shows the distribution of the Monte-Carlo particles within the computational domain. Assuring 
that this distribution is uniform allows us to discount any differences between LES-FD and VSFDF as 
resulting form interpolation bias. The influence of the average particle density, N E , on the first and second 
moments is shown in Figs. 4-12. It can be observed that N E does not significantly influence the first nor 
the second moments. The differences displayed on these figures seem visually comparable to those due 
to changing random number generator seed. The later are shown explicitly for different values of A E in 
Figs. 13-22. It is important to keep this influence in mind when dealing with Monte-Carlo type simulations 
as it is possible, but highly improbable, that two VSFDF simulations will yield identical results. Figs. 23- 
31 show the influence of the ensemble domain size (See Fig. 1) on the first and second moments. The 
first moments as obtained by VSFDF simulation agree very well with their LES-FD counterparts even for 
large values of A E ■ The small differences in these moments can be contributed to the difference in the 
seed of the random number generator for different cases. However, the second moments exhibit significant 
differences. These differences exhibit a converging pattern of behavior as the ensemble domain is decreased. 
This was predicted and shown theoretically in Sec. E. As ensemble domain size decreases, the values of 
the SGS components predicted by the VSFDF simulation converge to their LES-FD values. Figs. 16- 
19, and 25-28 illustrate this behavior quite clearly. The scatter plots of first and second moments are 
presented in Figs. 32-34. The correlation and linear regression coefficients (denoted by p and r respectively 
on the figures) offer another set of indicators for consistency for VSFDF simulations. The p 2 indicates the 
percentage of variance of the VSFDF data accounted for by a linear fit. Figs. 35-44 provide a simple visual 
demonstration of the consistency of the VSFDF simulations and Figs. 45 and 46 show the consistency of 
the time evolution of the simulation. 


4 Summary and Concluding Remarks 

The Filtered Density Function (FDF) methodology has proven very effective in LES of turbulent reactive 
flows. 24 ’ 25 ’ 27 ’ 28 ’ 30 In all of the previous investigations, however, the FDF of either only the scalar 24 (SFDF), 
or velocity 31 (VFDF) were considered. The objective of present work is to develop and validate the FDF 
methodology for LES of the joint velocity-scalar field (VSFDF). For this purpose, a methodology for 
velocity FDF developed by Gicquel 31 is extended to include the scalar field resulting in a VSFDF method 
for LES. The exact transport equations governing the evolution of VSFDF are derived. It is shown that 
the effects of the SGS convection in these equations as well as the terms resulting from reaction source 
terms appear in the closed form. The remaining unclosed terms are modeled via VSFDF formulation with 
Generalized Langevin Model (GLM) for the Lagrangian velocity field evolution and Linear Mean Squaie 
Estimation (LMSE) Model for the Lagrangian scalar field evolution. The closure strategy is similar to that 
for PDF methods in Reynolds Averaged Navier-Stokes (RANS) simulations and consequently equivalent 
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to a second-order moment SGS closure. 

The modeled VSFDF transport equations are solved numerically via a Lagrangian Monte Carlo scheme. 
The result of this simulation is then compared to a well established LES using a finite difference scheme 
(LES-FD). The consistency of the VSFDF method and the convergence of its Monte Carlo solutions are 
assessed. This assessment is done by comparing solutions of the first two filtered moments obtained via 
Monte Carlo procedure and LES-FD. The third order moment terms in the second filtered moment trans- 
port equations for LES-FD are closed by corresponding quantities obtained from VSFDF. The assessment 
itself contains the comparisons of Reynolds average quantities for varying ensemble domain size(A£) and 
average particle density siz q(Ne). The theoretical predictions for these comparisons are validated. Addi- 
tionally, instantaneous comparisons of data scatter, visual inspection of first and second moments, and time 
consistency of the simulations are examined. The consistency and convergence of the VSFDF simulation 
is demonstrated by good agreement of the first two moments as obtained with Monte Carlo simulations 
and LES-FD for all considered test cases. 

Work is in progress on further evaluating the performance of VSkDF model by comparing the results 
obtained with a 3D homogeneous shear layer case with the DNS and experimental data. The effectiveness of 
the model as compared to some commonly used LES models such as Smagorinsky, or Dynamic Smagor insky 
will also be evaluated as well its computational expense. Progress towards the development of V SFDF foi 
LES of reaction flows will follow as the VSFDFs major advantages over its predecessors are the treatment 
of the convective transport and reactive terms in a closed form. 
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Figure 4: Cross-stream variation of the Reynolds averaged values of the ( u ) component of the velocity 
field: (a) A £ = 0.5A, (b) A £ = 1.0, (c) A £ = 2.0 
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Figure 5: Cross-stream variation of the Reynolds averaged values of the (</>), scalar field: (a) A E - 0.5A 
(b) Ae = 1.0, (c) A e = 2.0 
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Figure 6: Cross-stream variation of the Reynolds averaged values of the t(u,u) component of the SGS 
stress tensor: (a) A^ = 0.5A, (b) A & = 1.0, (c) A# = 2.0 








Figure 7: Cross-stream variation of the Reynolds averaged values of the t{v,v) component of the SGS 
stress tensor: (a) Ag = 0.5A, (b) A e = 1-0, (c) A e — 2.0 
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Figure 8: Cross-stream variation of the Reynolds averaged values of the r(u,u), SGS stress tensor: (a) 
A E = 0.5A, (b) A E = 1-0, (c) A E = 2.0 
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Figure 9: Cross-stream variation of the Reynolds averaged values of the turbulent kinetic energy, k: (a) 
A E = 0.5A, (b) A E = 1.0, (c) A e = 2.0 







Figure 10: Cross-stream variation of the Reynolds averaged values of the t(u, 0), SGS tensor: (a) A e 
0.5A, (b) A E = 1.0, (c) A E = 2-0 
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Figure 11: Cross-stream variation of the Reynolds averaged values of the r(u,</>), SGS tensor: (a) A g 
0.5A, (b) A E = 1.0, (c) A e = 2.0 
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Figure 12: Cross-stream variation of the Reynolds averaged values of the SGS tensor: (a) A e 

0.5A, (b) A E = 1.0, (c) A e = 2.0 
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Figure 13: Effect of the random number sequence used in the Monte-Carlo solver with Ne — 40 01 
repeatability of the LES-FD and VSFDF simulations for Reynolds averaged values of (u) field. (— LES 
FD, VSFDF) 









0.01 
0 

- 0.01 

-20 -10 0 10 20 

Y 

Figure 14: Effect of the random number sequence used in the Monte-Carlo solver with N% = 40 on 
repeatability of the LES-FD and VSFDF simulations for Reynolds averaged values of (v) field. (- LES- 
FD, VSFDF) 
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Figure 15: Effect of the random nur 
repeatability of the LES-FD and VSF 
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Figure 16: Effect of the random number sequence used in the Monte-Carlo solver with Ne = 40 o 
repeatability of the LES-FD and VSFDF simulations for Reynolds averaged values of t(u,u ) field. ( 
LES-FD, VSFDF) 
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Figure 19: Effect of the random number sequence used in the Monte-Carlo solver with Ne = 40 on 
repeatability of the LES-FD and VSFDF simulations for Reynolds averaged values of turbulent kinetic 
energy, k , field. (— LES-FD, VSFDF) 









Figure 20: Effect of the random number sequence used in the Monte-Carlo solver with Ng — 40 c 
repeatability of the LES-FD and VSFDF simulations for Reynolds averaged values of t{u,4>) held. ( 
LES-FD, VSFDF) 










Figure 21: Effect of the random number sequence used in the Monte-Carlo solver with Ng — 40 o 
repeatability of the LES-FD and VSFDF simulations for Reynolds averaged values of t(v,4>) field. ( 
LES-FD, VSFDF) 








Figure 22: Effect of the random number sequence u; 
repeatability of the LES-FD and VSFDF simulations 
LES-FD, VSFDF) 
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Figure 23: Cross-stream variation of the Reynolds averaged values of the ( u ) component of the velocity 
field: (a) N E = 20, (b) N E = 40, (c) N E = 80 






Figure 24: Cross-stream variation of the Reynolds averaged values of the {</>), scalar field: (a) A'/.; 
(b) N e = 40, (c) N e = 80 
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Figure 25: Cross-stream variation of the Reynolds averaged values of the t(u,u ) component of the SGS 
stress tensor: (a) Ne = 20, (b) Ne = 40, (c) Ne = 80 
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Figure 26: Cross-stream variation of the Reynolds averaged values of the t(v,v ) component of the SGS 
stress tensor: (a) Ng = 20, (b) Ng = 40, (c) Ne = 80 
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Figure 27: Cross-stream variation of the Reynolds averaged values of the r(it, u), SGS stress tensor: (a) 
N e = 20, (b) N e = 40, (c) ' N e = 80 
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Figure 28: Cross-stream variation of the Reynolds averaged values of the turbulent kinetic energy, k : (a) 
N e = 20, (b) N e = 40, (c) N e = 80 
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Figure 29: Cross-stream variation of the Reynolds averaged values of the t(u,4>), SGS tensor: (a) iV# — 20 
(b) N e = 40, (c) N e = 80 
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Figure 30: Cross-stream variation of the Reynolds averaged values of the r(u, </>), SGS tensor: (a) N E 
(b) N e = 40, (c) N e = 80 
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Figure 31: Cross-stream variation of the Reynolds averaged values of the r((f>, <fi), SGS tensor: (a) Ne 
(b) Ne = 40, (c) Ne = 80 
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Figure 32: Scatter plot of the filtered velocity, (it) and ( v ), and scalar, (</>) obtained via VSFDF and 
LES-FD for A# = 0.5A,7V£; = 40. (a) ( u ), (b) ( v ), (c) (</)). - - least-squares-fit line, - 1:1 fit line, p is the 
correlation coefficient, r is the linear regression coefficient 
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Figure 33: Scatter plot of the components of the SGS stress tensor, r(u,u), t(v,v), and t(u,v) obtained 
via VSFDF and LES-FD for Ag = 0.5A ,N& = 40. (a) t(u,u), (b) r(v,v), (c) r(u,v). — least-squares-fit 
line, - 1:1 fit line, p is the correlation coefficient, r is the linear regression coefficient 
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Figure 34: Scatter plot of the components of the SGS tensor, r(«, 0 ), t(v, cj> ), and r(0, <f>) obtained via 
VSFDF and LES-FD for A# = 0.5A,Ne = 40. (a) t{u,4>), (b) t(u,<^>), (c) t(</>, </>). - - least-squares-fit 
line, - 1:1 fit line, p is the correlation coefficient, r is the linear regression coefficient 
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Figure 35: Visual consistency for stream-wise velocity field, (it). The LES-FD results are obtained for 
A e = 0.5A, N e = 40 at t=34.3. (Captions: Left - LES-FD, Right - VSFDF) 
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Figure 36: Visual consistency for span- wise velocity field, ( v ). The LES-FD results are obtained for 
A e = 0.5A, N e = 40 at t=34.3. (Captions: Left - LES-FD, Right - VSFDF) 
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Figure 39: Visual consistency for the r(u, u) component of the SGS stress tensor. The LES-FD results ar 
obtained for A E = 0.5A, N E = 40 at t=34.3. (Captions: Left - LES-FD, Right - VSFDF) 
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Figure 40: Visual consistency for the t(v,v) component of the SGS stress tensor. The LES-FD results are 
obtained for A E = 0.5A, N E = 40 at t=34.3. (Captions: Left - LES-FD, Right - VSFDF) 
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Figure 41: Visual consistency for the t(u,v) component of the SGS stress tensor. The LES-FD results 
obtained for Ag = 0.5A, N E = 40 at t=34.3. (Captions: Left - LES-FD, Right - VSFDF) 
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Figure 42: Visual consistency for the t(u, <f>) component of the SGS stress tensor. The LES-FD results 
obtained for A E = 0.5A, N E = 40 at t=34.3. (Captions: Left - LES-FD, Right - VSFDF) 
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Figure 43: Visual consistency for the t(v, <p) component of the SGS stress tensor. The LES-FD results 
obtained for A E = 0.5A, N E = 40 at t=34.3. (Captions: Left - LES-FD, Right - VSFDF) 



Figure 44: Visual consistency for the r(<p, </>) component of the SGS stress tensor. The LES-FD results 
obtained for A E = 0.5A, N E = 40 at t=34.3. (Captions: Left - LES-FD, Right - VSFDF) 
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Figure 45: Time evolution of the scalar(with superimposed vorticity isolines) and vorticity fields for LES- 
FD obtained with A # = 0.5A and Ng = 40 
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Figure 46: Time evolution of the 
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lposed vorticity isolines) and vorticity fields for VSFDF 
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